/*
 * Academic License - for use in teaching, academic research, and meeting
 * course requirements at degree granting institutions only.  Not for
 * government, commercial, or other organizational use.
 *
 * linearLeastSquares.c
 *
 * Code generation for function 'linearLeastSquares'
 *
 */

/* Include files */
#include "linearLeastSquares.h"
#include "rt_nonfinite.h"
#include "xnrm2.h"
#include "xzlarf.h"
#include "xzlarfg.h"
#include <math.h>

/* Function Definitions */
void linearLeastSquares(double lhs[32], double rhs[8], double dx[4])
{
  double jpvt[4];
  double work[4];
  double d;
  double temp;
  int b_i;
  int i;
  int ii;
  int ix;
  int iy;
  int j;
  int k;
  int nfxd;
  int temp_tmp;
  jpvt[0] = 0.0;
  jpvt[1] = 0.0;
  jpvt[2] = 0.0;
  jpvt[3] = 0.0;
  nfxd = 0;
  for (j = 0; j < 4; j++) {
    if (jpvt[j] != 0.0) {
      nfxd++;
      if (j + 1 != nfxd) {
        ix = j << 3;
        iy = (nfxd - 1) << 3;
        for (k = 0; k < 8; k++) {
          temp_tmp = ix + k;
          temp = lhs[temp_tmp];
          i = iy + k;
          lhs[temp_tmp] = lhs[i];
          lhs[i] = temp;
        }
        jpvt[j] = jpvt[nfxd - 1];
        jpvt[nfxd - 1] = (double)j + 1.0;
      } else {
        jpvt[j] = (double)j + 1.0;
      }
    } else {
      jpvt[j] = (double)j + 1.0;
    }
    dx[j] = 0.0;
    work[j] = 0.0;
  }
  i = (unsigned char)nfxd;
  for (b_i = 0; b_i < i; b_i++) {
    ii = (b_i << 3) + b_i;
    temp = lhs[ii];
    d = xzlarfg(8 - b_i, &temp, lhs, ii + 2);
    dx[b_i] = d;
    lhs[ii] = temp;
    if (b_i + 1 < 4) {
      lhs[ii] = 1.0;
      xzlarf(8 - b_i, 3 - b_i, ii + 1, d, lhs, ii + 9, work);
      lhs[ii] = temp;
    }
  }
  if (nfxd < 4) {
    double vn1[4];
    double vn2[4];
    work[0] = 0.0;
    vn1[0] = 0.0;
    vn2[0] = 0.0;
    work[1] = 0.0;
    vn1[1] = 0.0;
    vn2[1] = 0.0;
    work[2] = 0.0;
    vn1[2] = 0.0;
    vn2[2] = 0.0;
    work[3] = 0.0;
    vn1[3] = 0.0;
    vn2[3] = 0.0;
    i = nfxd + 1;
    for (j = i; j < 5; j++) {
      d = xnrm2(8 - nfxd, lhs, (nfxd + ((j - 1) << 3)) + 1);
      vn1[j - 1] = d;
      vn2[j - 1] = d;
    }
    for (b_i = i; b_i < 5; b_i++) {
      double s;
      int ip1;
      ip1 = b_i + 1;
      j = (b_i - 1) << 3;
      ii = (j + b_i) - 1;
      iy = 5 - b_i;
      nfxd = -1;
      if (5 - b_i > 1) {
        temp = fabs(vn1[b_i - 1]);
        for (k = 2; k <= iy; k++) {
          s = fabs(vn1[(b_i + k) - 2]);
          if (s > temp) {
            nfxd = k - 2;
            temp = s;
          }
        }
      }
      iy = b_i + nfxd;
      if (iy + 1 != b_i) {
        ix = iy << 3;
        for (k = 0; k < 8; k++) {
          temp_tmp = ix + k;
          temp = lhs[temp_tmp];
          nfxd = j + k;
          lhs[temp_tmp] = lhs[nfxd];
          lhs[nfxd] = temp;
        }
        temp = jpvt[iy];
        jpvt[iy] = jpvt[b_i - 1];
        jpvt[b_i - 1] = temp;
        vn1[iy] = vn1[b_i - 1];
        vn2[iy] = vn2[b_i - 1];
      }
      temp = lhs[ii];
      d = xzlarfg(9 - b_i, &temp, lhs, ii + 2);
      dx[b_i - 1] = d;
      lhs[ii] = temp;
      if (b_i < 4) {
        lhs[ii] = 1.0;
        xzlarf(9 - b_i, 4 - b_i, ii + 1, d, lhs, ii + 9, work);
        lhs[ii] = temp;
      }
      for (j = ip1; j < 5; j++) {
        iy = b_i + ((j - 1) << 3);
        d = vn1[j - 1];
        if (d != 0.0) {
          temp = fabs(lhs[iy - 1]) / d;
          temp = 1.0 - temp * temp;
          if (temp < 0.0) {
            temp = 0.0;
          }
          s = d / vn2[j - 1];
          s = temp * (s * s);
          if (s <= 1.4901161193847656E-8) {
            d = xnrm2(8 - b_i, lhs, iy + 1);
            vn1[j - 1] = d;
            vn2[j - 1] = d;
          } else {
            vn1[j - 1] = d * sqrt(temp);
          }
        }
      }
    }
  }
  for (j = 0; j < 4; j++) {
    if (dx[j] != 0.0) {
      temp = rhs[j];
      i = j + 2;
      for (b_i = i; b_i < 9; b_i++) {
        temp += lhs[(b_i + (j << 3)) - 1] * rhs[b_i - 1];
      }
      temp *= dx[j];
      if (temp != 0.0) {
        rhs[j] -= temp;
        for (b_i = i; b_i < 9; b_i++) {
          rhs[b_i - 1] -= lhs[(b_i + (j << 3)) - 1] * temp;
        }
      }
    }
  }
  dx[0] = rhs[0];
  dx[1] = rhs[1];
  dx[2] = rhs[2];
  dx[3] = rhs[3];
  for (j = 3; j >= 0; j--) {
    iy = j + (j << 3);
    dx[j] /= lhs[iy];
    for (b_i = 0; b_i < j; b_i++) {
      ix = (j - b_i) - 1;
      dx[ix] -= dx[j] * lhs[(iy - b_i) - 1];
    }
  }
  work[1] = dx[1];
  work[2] = dx[2];
  work[3] = dx[3];
  dx[(int)jpvt[0] - 1] = dx[0];
  dx[(int)jpvt[1] - 1] = work[1];
  dx[(int)jpvt[2] - 1] = work[2];
  dx[(int)jpvt[3] - 1] = work[3];
}

/* End of code generation (linearLeastSquares.c) */
